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Abstract 


We  describe  an  algorithm  for  the  solution  of  the  laminar,  incompressible  Navier- 
Stokes  equations.  The  basic  algorithm  is  a  multigrid  method  based  on  a  robust, 
box-based  smoothing  step.  Its  most  important  feature  is  the  incorporation  of 
automatic,  dynamic  mesh  refinement.  Using  an  approximation  to  the  local  trun¬ 
cation  error  to  control  the  refinement,  we  use  a  form  of  domain  decomposition  to 
introduce  patches  of  finer  grid  wherever  they  arc  needed  to  ensure  an  accurate 
solution.  This  refinement  strategy  is  completely  local:  regions  that  satisfy  our 
tolerance  arc  unmodified,  except  when  they  must  be  refined  to  maintain  reason¬ 
able  mesh  ratios.  This  locality  has  the  important  consequence  that  boundary 
layers  and  other  regions  of  sharp  transition  do  not  “steal”  mesh  points  from  sur¬ 
rounding  regions  of  smooth  flow,  in  contrast  to  moving  mesh  strategies  where 
such  “stealing”  is  inevitable. 


Our  algorithm  supports  generalized  simple  domains,  that  is,  any  domain  defined 
by  horizontal  and  vertical  lines.  This  generality  is  a  natural  consequence  of  our 
domain  decomposition  approach.  We  base  our  program  on  a  standard  staggered- 
grid  formulation  of  the  Navier-Stokes  equations  for  robustness  and  efficiency. 
To  ensure  discrete  mass  conservation,  we  have  introduced  special  grid  transfer 
operators  at  grid  interfaces  in  the  multigrid  algorithm.  While  these  operators 
complicate  the  algorithm  somewhat,  our  approach  results  in  exact  mass  conserva¬ 
tion  and  rapid  convergence. 


In  this  paper,  we  present  results  for  three  model  problems:  the  driven-cavity,  a 
backward-facing  step,  and  a  sudden  expansion/contraction.  Our  approach  obtains 
convergence  rates  that  are  independent  of  grid  size  and  that  compare  favorably 
with  other  multigrid  algorithms.  We  aLso  compare  the  accuracy  of  our  results 
with  other  benchmark  results. 
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1.  Introduction 

Computational  fluid  dynamics  is  one  of  several  fields  where  one  finds  boundary  and  interior 
layers  in  complex  interactions.  These  regions  of  sharp  transition  frequently  play  crucial  roles  in 
determining  the  global  solution,  despite  their  small  sizes.  The  need  for  true  local  refinement  is 
therefore  great. 

We  use  the  phrase  “true  local  refinement”  to  mean  the  introduction  of  extra  unknowns  in 
areas  of  interest  and  only  in  tliose  areas.  This  contrasts  with  approaches  that  move  mesh  points 
into  areas  of  interest  but  maintain  a  constant  number  of  unknowns.  Our  strategy  is  atypical  of 
finite  difference-based  codes,  including  those  that  u.se  mesh  transformations,  since  arrays  arc  gen¬ 
erally  the  only  data  structure  used.  While  moving  mesh  points  can  improve  accuracy,  this 
approach  gives  limited  control  over  global  accuracy.  The  necessary  flexibility  is  easily  provided 
by  finite  clement  codes,  but  the  overhead  in  mt,’ntaining  this  freedom  slows  solution  speed.  The 
essential  issue  is  one  of  grid  structure  and  regularity.  We  enhance  a  regular  mesh  with  the 
pointers  necessary  u  allow  us  to  refine  specific  sections  of  the  grid  while  leaving  the  remainder 
unmodified. 

Our  approximat'on  is  based  on  the  use  of  finite  differences  at  each  grid  level.  We  extend  the 
usual  formulation  to  support  locally  refined  patches.  Such  an  extension  requires  a  data  structure 
that  is  simple  enough  to  implement  efficiently,  yet  flexible  enough  to  permit  the  very  complicated 
patterns  of  local  mesh  refinement  that  would  be  generated  by  an  adaptive  mesh  refinement 
strategy.  To  achieve  these  goals,  we  have  selected  a  quad-tree  data  structure.  At  any  level,  the 
grid  is  generated  by  a  collection  of  square  grid  patches  containing  cells,  where  N  =  2p  and  p 
is  a  specified  integer  greater  than  one.  The  collection  of  patches  making  up  a  grid  may  or  may 
not  be  contiguous.  A  grid  refinement  is  generated  by  selectively  refining  any  or  all  of  the  qua¬ 
drants  in  any  of  the  patches  in  tlie  grid.  Refinement  by  quadrants  allows  finer-grained  adaptivity 
than  refining  only  entire  patches.  The  refined  quadrant  becomes  a  patch  at  the  next  finer  level, 
with  the  same  structure  as  its  parent  patch.  Thus  implementation  is  simple,  yet  the  scheme  is 
flexible  enough  to  allow  coupling  with  adaptive  mesh  slialegies.  In  addition,  this  quad-tree  data 
structure  allows  the  mulligrid  algorithm  to  be  easily  adapted  to  shared-memory  parallel  architec¬ 
tures  [15). 

Any  multigrid  algorithm  based  on  the  use  of  nonuniform  composite  grids  must  address  two 
numerical  is.sucs.  First  is  the  issue  of  approximating  the  governing  equations  on  a  nonuniform 
composite  grid.  In  particular,  we  seek  consistent  approximations  that  do  not  violate  the  conserva¬ 
tion  laws  beyond  the  intrinsic  level  of  the  approximation.  Second,  the  nature  of  the  transfer  func¬ 
tions  on  a  composite  grid  must  be  specified  so  that  a  consistent  multigrid  approximation  results. 

Our  algorithms  were  tested  on  model  problems  in  two  dimensions;  in  this  paper  we  pre.sent 
results  on  a  driven-cavity,  the  backward-facing  step,  and  a  sudden  expansion/contraction 
configuration.  These  problems  are  sufficiently  demanding  to  test  the  numerical  procedures,  while 
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being  well  understood  with  an  abundance  of  benchmark  data.  The  point  is  to  design  algorithms 
that  easily  extend  to  more  complex  geometries  and  more  dilTicult  flows. 

The  essential  features  of  our  approach  are  a  dynamically  adaptive  grid  composed  of 
dilTercnt-sizcd  rectangles  and  a  discretization  based  on  hybrid  differences,  with  an  interface  treat¬ 
ment  that  conserves  mass  exactly.  The  solution  algorithm  uses  a  multigrid  approach  using 
“boxed-based”  relaxation  and  a  variety  of  cycling  strategies.  This  is  similar  to  the  method 
adopted  by  Thompson  and  Ferziger  [15].  However,  there  are  a  number  of  important  differences. 
Our  adoption  of  “patches”  as  our  basic  data  object  and  our  di-i^erent  choice  of  grid  partitioning 
have  allowed  us  to  avoid  the  step  of  grouping  cells  marked  for  refinement  into  subrectangles 
(which  is  potentially  costly  [15]).  Moreover,  we  can  handle  somewhat  more  general  geometries, 
support  more  complex  refinement  patterns,  and  (as  we  show  in  reference  [16])  efficiently  perform 
calculations  in  parallel. 

Our  experiments  have  shown  that  the  novel,  mass-conserving  transfer  operators  (which  arc 
also  an  integral  part  of  our  interface  treatment)  are  both  natural  to  use  and  accurate.  The  use  of 
xlf  as  an  error  estimator  works  well  and  is  also  natural  in  the  multigrid  context. 

2.  The  Test  Problems  and  Governing  Equations 

We  consider  steady,  two-dimensional,  laminar,  incompressible  flow.  In  this  section  we 
describe  our  simplest  test  problem,  namely,  flow  in  a  square  box  that  is  driven  by  a  moving  lid 
(the  so-called  driven-cavity  problem).  In  Section  8  we  describe  the  other  two  geometries  we  have 
considered  (the  backward-facing  step  and  the  sudden  expansion/contraction)  in  the  context  of  our 
results.  The  code  can  easily  handle  any  domain  defined  by  horizontal  and  vertical  lines  by  a 
natural  extension  of  our  technique  for  handling  grid  refinement.  The  problems  are  well  known; 
we  have  chosen  them  bccau.se  they  are  geometrically  simple  but  display  complex  flow  stmetures. 
Also,  there  arc  many  well-documented  numerical  solutions  (see,  for  example,  [5],  [2],  [7],  and 
[I?]). 

Our  governing  equations  arc 

where  the  Reynolds  number  is  defined  by 
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Re  =  ^  L/v. 

We  take  the  length  scale  to  be  L  =  1;  w  is  the  characteristic  velocity  to  be  normalized  to  one.  For 
the  driven-cavity  problem,  the  domain  is  the  unit  square,  and  the  associated  boundary  conditions 
are 

«  =  0  ,  V  =  0  on  the  bottom  and  sides  (y  =  0,  x  =0  and  ;«:  =  1) 
and 

u  =  \  ,  V  =  0  on  the  top  (y  =  1) 

(see  Figure  1(a)).  Other  problem  domains  and  boundary  conditions  are  shown  in  Section  8. 

Equations  (1)  represent  a  set  of  three  coupled  nonlinear  partial  differential  equations.  The 
degree  of  nonlinearity  is  determined  oy  the  Reynolds  number  (Re).  We  present  two  sets  of 
results;  one  at  Re  =  I,  which  is  almost  linear  but  allows  us  to  lor>k  at  certain  asymptotic  properties 
of  the  discretization,  and  one  at  Re  =  400,  which  is  a  .casonably  nonlinear  problem  but  stili 
corresponds  to  stable,  steady,  physical  solutions.  Of  course,  for  this  test  problem,  other  formula¬ 
tions,  such  as  stream-function  vorticity,  are  available.  We  have  retained  the  primitive  variable 
formulation  because  it  is  easier  to  generalize  to  more  complex  flow  regimes  and  to  three- 
dimensional  problems. 

In  addition  to  the  two  comer  singularities  mentioned,  there  is  the  usual  pressure  indeter¬ 
minacy,  always  found  with  incompressible  flows:  the  pressure  is  determined  only  up  to  an  addi¬ 
tive  constant.  Since  our  algorithm  is  ba.sed  on  iterative  methods,  we  simply  ignore  the  additive 
constant.  Our  iterative  scheme  converges  to  a  solution  with  average  pressure  almost  the  same  as 
that  of  the  starting  guess.  We  avoid  some  unnecessary  complexity,  and  we  gain  in  convergence 
rate  by  this  approach. 


3.  Basic  Discretization 

We  have  chosen  a  standard  staggered  grid  as  the  basic  discretization  scheme.  In  this  formu¬ 
lation  the  computational  domain  is  filled  with  “mass  control  volumes.”  Pressures  are  defined  in 
the  centers  of  these  rectangles,  vertical  velocity  components  are  defined  at  the  top  and  bottom 
faces,  and  horizontal  ones  arc  similarly  offset  (see  Figure  1(b)).  This  means  that  it  is  not  neces¬ 
sary  to  generate  pressure  boundary  conditions. 

We  have  several  reasons  for  choosing  this  discretization.  Principally,  the  properties  of  this 
approach  arc  well  understood,  and  we  wish  to  avoid  possible  interactions  between  newer  discreti¬ 
zations  and  our  treatment  of  the  finc/coansc  interfaces.  The  main  advantage  is  the  stability  of  the 
scheme:  it  will  not  exhibit  spurious  pre.ssurc  modes,  provided  only  that  the  momentum  equations 
are  stably  discretized.  Moreover,  the  discrete  cilipticity  is  somewhat  better  than  in  nonstaggered 
formulations,  although  these  can  be  corrected  by  introducing  artificial  compressibility  into  the 
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continuiry  cqualicn. 

Our  approximations  of  derivatives  are  also  quite  standard.  We  treat  diffusion  terms  with 
standard  three-point  central  ditTercnces.  The  pressure  gradient  terms  are  also  discretized  by  using 
central  differences.  The  convection  terms  require  more  careful  handling.  We  use  upwind 
differencing  if  the  mesh  Reynolds  number  (uhtv)  is  greater  than  2;  otherwise,  we  use  central 
differencing.  Thus,  for  sufficiently  small  meshes,  this  discretization  has  second-order  local  trun¬ 
cation  errors  for  all  variables.  An  equally  important  property  is  that  the  discretization  is  stable 
(/i -elliptic)  for  all  flow  rates  on  all  grids  [13|.  This  is  essential  in  multigrid  methods:  violation  of 
this  condition  can  mean  that  the  coarsc-grid  correction  is  grossly  inaccurate.  The  actua'  discrete 
approximations  are  given  in  detail  in  [1  j. 

Mesh  transformations  have  become  an  increasingly  common  technique  for  adapting  finite 
difference  algorithms  to  complex  geometries.  Nonstaggered  discretizations  have  an  obvious 
advantage  in  this  situation,  since  it  is  necessary  to  calculate  and  store  only  one  set  of  geometric 
information  instead  of  fo  ir  (in  three  dimensions).  However,  the  amount  of  increased  computa¬ 
tional  complexity  in  writing  the  differential  equations  in  transformed  coordinates  is  very  large, 
and  it  is  not  clear  whether  this  is  the  best  method  for  dealing  with  complex  geometries.  One 
approach  that  fits  naturally  into  our  multigrid  framework  is  to  have  locally  fitted  grids  that  over¬ 
lap  with  a  Cartesian  system  in  the  interior  (see,  for  example  Chesshire  and  Henshaw  [6]).  This 
method  allows  efficient  treatment  of  the  interiors  while  accurately  representing  some  classes  of 
complex  geometry. 


4.  A  Brief  Overview  of  the  Multigrid  Solution  Algorithm 

In  the  current  context,  the  goal  of  a  good  multigrid  algorithm  can  be  defined  as  the  robust 
solution  of  a  set  of  elliptic  equations  in  an  amount  of  work  equivalent  to  a  small  number  of  relax¬ 
ation  sweeps.  In  particular,  we  reject  strategies  where  the  convergence  rates  degrade  as  the 
meshes  become  finer.  Section  8  describes  the  specific  aims  of  our  project  in  more  detail. 

For  completeness,  we  give  a  brief  outline  of  our  multigrid  algorithm  here.  A  more  com¬ 
plete  description  of  multigrid  methods  is  given  in  references  [4],  [9],  and  [10];  and  a  description 
of  an  algorithm  similar  to  the  present  one,  but  for  nonadapted  grids,  is  given  in  [1].  Thus  we  res¬ 
trict  ourselves  here  to  an  overview  of  our  current  algorithm,  stressing  new  features  where 
appropriate. 

We  have  already  defined  the  fine-grid  representation  of  our  problem.  For  efficiency,  we  use 
geometric  coarse-grid  operators.  (In  other  words,  the  coarse-grid  operator  is  defined  in  the  same 
way  as  the  fine-grid  one  except  that  h.  is  replaced  by  2h .) 

The  staggered  di.scrctization  works  well  in  the  context  of  multigrid  methods  and  offers  a 
number  of  benefits.  Among  the  various  ways  of  performing  the  relaxations,  we  use  the  approach 
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of  simulianeously  updating  all  five  variables  (four  velocity  components  and  pressure)  associated 
with  a  mass  control  volume.  This  method  seems  to  be  more  effective  and  more  reliable  than 
schemes  that  cyclically  update  the  variables  associated  with  each  cell,  especially  at  high  Reynolds 
numbers.  At  fine/coarse-grid  interfaces,  the  staggered  formulation  has  the  advantage  that  the 
molecules  that  must  be  simultaneously  updated  are  smaller  than  those  in  nonstaggered  formula¬ 
tions;  hence,  accuracy  may  be  improved.  Details  are  given  later. 

Before  the  update  of  each  cell,  we  locaUy  linearize  about  the  old  solution.  There  are  some 
indications  that  we  would  achieve  better  rates  of  convergence  if  we  used  local  Newton  updates. 
However,  the  increased  amount  of  arithmetic  involved  in  this  process  offsets  the  faster  conver¬ 
gence.  On  the  coarsest  grid  we  perform  as  many  sweeps  as  are  necessary  to  satisfy  a  convergence 
criterion.  On  finer  grids  we  perform  a  small  number  of  sweeps,  as  defined  by  our  cycling  strategy. 

We  use  full  weighting  for  the  multigrid  restrictions  in  all  cases.  This  gives  us  better  conver¬ 
gence  properties,  since  it  minimizes  aliasing  effects.  (This  technique  contrasts  with  our  earlier 
approach  [1],  where  we  used  linear  interpolation  for  our  restrictions.) 

For  prolongations  of  velocities  we  use  an  approach  thai  is  second  order  but  automatically 
satisfies  the  discrete  continuity  equations  on  fine  cells,  provided  they  were  satisfied  on  the 
corresponding  coarse-grid  ccU.  This  property  somewhat  simplifies  our  treatment  of  interfaces. 
Details  are  given  in  the  next  section.  Pre.s.sitres  arc  prolonged  by  using  linear  interpolation.  It 
should  be  noted  that  the  proof  given  in  [1],  showing  that  discrete  continuity  on  the  coarsest  grid 
implies  discrete  continuity  on  aU  grids,  is  still  valid  in  this  context. 

We  have  implemented  all  the  common  multigrid  cycling  strategies  (F,  W,  and  V)  as  well  as 
adaptive  cycling.  Earlier  results  have  shown  that  for  nonlinear  Navier-Stokes  equations,  F  cycles 
are  somewhat  more  efficient  [1].  We  have  done  most  of  our  work  with  this  strategy.  However, 
experiments  have  shown  that  the  rate  of  convergence  is  insensitive  to  the  cycling  strategy.  Note 
that  our  cycling  strategy  is  the  same  in  the  uniform  and  nonuniform  grid  case;  we  perform  opera¬ 
tions  on  cells  of  a  particular  mesh  size  during  each  sweep. 


5.  The  Discretization  on  Compound  Grids  (Allowable  Grids) 

It  is  common,  especially  in  multigrid  programs,  to  use  uniform  grids  at  each  level  of 
refinement.  Such  programs  are  easier  to  write  and  avoid  the  computational  overhead  of  dealing 
with  varying  mesh  spacing.  Of  course,  in  many  problems  uniform  grids  lead  to  excessive  prob¬ 
lem  size  in  order  to  achieve  a  desired  degree  of  resolution  in  the  solution.  Our  approach  is  to 
recover  the  flexibility  of  nonunifonn  grids  by  generating  a  composite  grid  made  up  of  a  collection 
of  patches  with  each  patch  having  a  uniform  grid.  The  size  of  each  patch  varies  according  to  its 
level  and  thereby  generates  a  nonuniform  composite  grid. 
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To  illustrate  this  concept,  let  us  start  with  patches  composed  of  4  cells.  (A  patch  consists  of 
A 2  cells,  where  N  =  2p ;  in  our  code  we  assume  that  p  >  1,  but  we  use  p  =  1  in  the  diagram  for 
clarity.)  Consider  a  square  domain  (such  as  a  driven  cavity)  with  a  composite  grid  made  up  of 
three  grid  levels.  We  shall  generate  a  typical  composite  grid  to  illustrate  the  procedure.  Suppose 
grid  level  1  consists  of  one  patch  that  is  the  entire  domain,  as  shown  in  Figure  2(a).  In  this  figure 
we  have  indicated  the  positions  associated  with  each  of  the  6  horizontal  velocities,  6  vertical  velo¬ 
cities,  and  4  pressures  defined  on  a  staggered  grid  associated  with  a  4-cell  patch.  Now,  each  patch 
is  made  up  of  4  quadrants.  As  patches  for  grid  level  2,  we  can  select  any  one,  two,  three,  or  all  of 
these  quadrants  as  patches  on  level  2.  For  illustrative  purposes,  we  select  two  patches  on  grid 
level  2.  Let  these  patches  be  the  lower  left  (patch  2)  and  upper  right  (patch  3)  quadrant  of  the  sin¬ 
gle  patch  (patch  1)  on  level  1  as  shown  in  Figure  2(b).  Now  consider  the  third  grid  level.  As  can¬ 
didates  for  patches  on  this  level,  we  can  select  any  combination  of  quadrants  in  patches  2  and  3. 
For  the  sake  of  simplicity  in  the  figures,  suppose  we  select  one  patch  at  level  3  located  in  the 
upper  right  quadrant  of  patch  3.  This  is  the  fourth  patch  in  the  composite  grid.  We  note  that  if 
we  had  selected,  for  example,  the  lower  left  quadrant  of  patch  3  to  be  a  patch  on  level  3,  we 
would  have  added  four  more  patches  on  level  2  so  that  adjacent  patches  are  separated  by  at  most 
one  level  (this  includes  diagonal  neighbors).  In  Figure  2  we  show  the  composite  grid  at  each 
level,  with  the  velocity  and  pressure  node  locations  indicated.  Thus,  the  final  composite  3  level 
grid  is  made  up  of  4  patches  containing  18  horizontal  velocity,  18  vertical  velocity,  and  13  pres¬ 
sure  nodes.  Note  that,  for  illustrative  purposes,  we  have  used  4-cell  patches;  in  practice  we  use 
16-ccll  (or  larger)  patches. 

Having  illustrated  the  generation  of  nonuniform  composite  grids  by  means  of  patches  gen¬ 
erated  by  quadrant  refinement,  we  now  consider  the  issues  involved  in  approximating  the  Navier- 
Stokes  equation  on  nonuniform  grids.  In  Figure  3(a)  we  show  a  typical  interface  between  grids  of 
different  sizes.  For  example.  Figure  3(a)  could  illustrate  an  interface  in  Figure  2(b)  between  a 
level- 1  grid  and  a  level-2  grid,  or  it  could  illustrate  an  interface  in  Figure  2(c)  between  a  level-2 
grid  and  a  levcl-3  grid. 

There  are  two  coarse-grid  cells.  The  right-hand  one  is  subdivided  into  four  fine-grid  cells. 
The  staggering  of  the  variables  is  shown  in  this  diagram.  We  have  shown  the  locations  of  all 
genuine  variables  in  Figure  3(a).  (By  “genuine”  we  mean  unknowns  that  are  associated  with  an 
approximation  to  some  differential  equation.)  We  write  a  horizontal  momentum  equation  for 
each  u  -velocity  (the  vertical  component  is  treated  similarly);  continuity  equations  are  written  at 
the  centers  of  cells  where  the  pressure  is  defined. 

The  usual  five-point  molecules  that  are  used  away  from  interfaces  must  be  modified  here. 
For  example,  the  continuity  equation  for  the  left-hand  cell  involves  three  « -components;  the 
momentum  equations  arc  also  modified  to  involve  both  coarse-  and  fine-grid  quantities.  Since  this 
approach  is  rather  inconvenient,  we  introduce  a  border  of  virtual  cells  in  the  coarse-grid  region 
(sec  F'g'T''  ^(3)).  These  arc  of  si^f’  ^  (ihe  fine-grid  spacing).  In  this  ca«e  seven  ^xtra  unknowns 
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are  associated  with  the  fine  grid.  Additionally,  four  extra  unknowns  are  associated  with  the 
coarse  grid.  Both  these  sets  of  unknowns  are  introduced  as  a  computational  convenience  to  “reg¬ 
ularize”  the  structure  of  unknowns.  Effectively,  the  new  fine-grid  cells  act  as  boundary  values  for 
the  fine-grid  equations,  and  those  for  the  coarse  grid  have  a  similar  role.  However,  we  have  to 
introduce  new  matching  conditions  that  are  consistent  with  the  differential  equations  but  that  are 
not  necessarily  derived  from  them  directly.  Obviously,  some  extra  work  is  involved  in  this  pro¬ 
cess,  but  in  our  multigrid  context  this  work  is  small  and  can  be  easily  incorporated  into  our  nor¬ 
mal  transfer  operators. 

For  efficiency  reasons  we  choose  to  relax  only  those  cells  that  have  a  particular  mesh  size 
(i.e.,  at  a  given  grid  level)  at  any  one  sweep.  We  relax  over  the  whole  domain  only  on  those 
coarse  grids  that  cover  the  whole  region.  Where  we  specify  refined  subdomains,  we  apply  the 
smoother  only  to  that  subgrid.  This  approach  follows  the  multilevel  adaptive  technique  of  Bai 
and  Brandt  [11]  and  the  fully  adaptive  composite  approach  of  McCormick  [12). 

We  have  constructed  this  arrangement  to  allow  a  “natural”  treatment  of  all  the  genuine  vari¬ 
ables.  This  means  that  we  write  the  usual  discretized  partial  differential  equations  for  all  the  vari¬ 
ables  shown  in  Figure  3(a).  All  of  these  variables  are  updated  by  our  usual  relaxation  technique 
throughout  the  smoothing  phase  of  the  multigrid  cycle  for  the  appropriate  grid  level  (see  jl]  for 
the  actual  formulae  and  relaxation  algorithm). 

The  dummy  unknowns  (shown  in  Figure  3(b))  are  not  updated  by  the  smoothing  operation. 
These  values  are  determined  by  interpolation  when  the  variables  are  transferred  from  the  finer  or 
coarser  grid  (dctx;nding  upon  the  stage  of  tbf  mi'itigrid  cvclingV  These  values  are  set  at  the 
beginning  of  the  relaxation  phase  and  act  as  boundary  values  while  the  genuine  variables,  which 
are  now  separated  from  the  interface,  are  updated.  The  interpolation  formulae  used  to  define  these 
dummy  variables  are  given  in  the  next  section.  Note,  however,  that  these  latter  variables  are 
updated  m  the  course  of  our  multigrid  cycling  so  that,  at  convergence,  both  our  discretized 
differential  equations  and  our  interface  conditions  arc  satisfied  simultaneously. 


6.  Interface  Conditions 

Our  standard  transfer  operators  use  linear  interpolation  for  all  variables.  This  works  reason¬ 
ably  well,  although  for  optimal  performance  we  should  use  higher-order  interpolation  the  first 
time  we  visit  a  new  grid  [9].  However,  our  numerical  experiments  show  that  our  strategy  works 
well  when  there  is  no  grid  refinement,  and  it  is  reasonably  simple  to  implement. 

If  there  is  no  grid  refinement,  the  transfer  operators  do  not  affect  Uie  accuracy  of  the  discreti¬ 
zation;  at  convergence  wc  solve  the  nonlinear,  fine-grid  equations  independent  of  the  transfers.  At 
interfaces,  however,  the  dummy  variables  do  affect  the  solution  because  they  are,  implicitly, 
boundary  values  for  the  two  subdomains.  It  is  natural  to  generate  these  values  bv  using  standard 
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transfer  operators.  However,  this  approach  can  lead  to  two  types  of  difficulty. 

The  primary  consideration  is  that  there  is  no  solution  to  the  di.scrcte  equatioas  unless 
discrete  mass  is  conserved.  If  we  consider  the  Hux  across  the  (horizontal)  interface,  we  require 

lh\'(j)  =  hvij  - \)  +  hv(j-  +  1). 

where  we  have  assumed  tliat  the  fine  grid  size  is  h  and  that  the  nomenclature  is  as  defined  in  Fig¬ 
ure  4.  This  requirement  is  automatically  satisfied  by  our  standard  re.striclion  operation.  However, 
nonnal  linear  prolongation  does  not  have  this  property.  We  have  modified  our  prolongation 
operator  .so  that  mass  is  conserved  by  this  operator  as  well.  We  define  the  line-grid  neighbors  of  a 
coarse-grid  velocity  by 

V  {j  +  ])  -  V{j )  -t-  {V  (J  -t-  1 )  —  V^(j  —  1  ))/8, 
and 

v(J  -  1)  =  +  \)-V(J  -  l))/8. 

These  expressions  are  formally  second-order  accurate,  and  insure  satisfaction  of  the  continuity 
equation  across  the  interface.  It  is  intere.sting  to  note  that  there  are  no  interpolation  formulae 
based  on  Lagrange  polynomials  that  satisfy  the  discrete  con.scrvation  property  and  that  have  order 
higher  than  two. 

The  other  aspect  of  the  interface  treatment  is  its  effect  on  accuracy.  Talcing  finite  differences 
of  values  obtained  by  interpolation  is,  in  general,  rather  suspect.  This  is  exactly  what  we  do, 
since  the  momentum  equation  at  the  interface  involves  differences  of  “dummy  variables” 
obtained  by  intcipolation.  One  can  show  that  with  the  linear  interpolation  our  treatment  of  advcc- 
uoii  tcmac  is  first-order  accurate  at  interfaces,  while  the  treatment  of  viscosity  terms  is  zeroth- 
order  accurate.  Global  truncation  order  is  ordinarily  one  order  of  accuracy  belter,  implying  global 
second-order  accuracy  for  advcction-dominalcd  problems,  but  only  first-order  accuracy  when 
viscosity  terms  dom.inate. 

Quadratic  interpolation  improves  accuracy  by  one  power  of  h,  giving  consistent  treatment 
of  viscosity  terms  as  well  at  interfaces.  Thus,  "'ilh  quadratic  interpolation,  we  can  expect 
.second-order  global  accuracy,  even  when  vi.scosity  terms  dominate,  a.s.suming  there  are  only  a 
bounded  number  of  interfaces  in  the  domain.  Even  in  the  case  where  the  number  of  interfaces 
grow's  w'ithout  bound,  advection-dominated  problems  remain  second-order  accurate.  All  of  this 
.suggests  that  the  quadratic  interpolation  is  .superior;  in  practice,  one  sees  relatively  little 
difference.  The  reason  seems  to  be  that  for  the  advection-dominated  problems  considered,  treat¬ 
ment  of  viscosity  terms  is  relative'y  unimportant.  The  simple  linear  interpolation  scheme  already 
gives  global  second-order  accuracy  in  most  cases  of  interest  and  was  used  in  most  of  our  experi¬ 
ments,  since  it  is  computationally  more  efficient.  However,  it  is  important  to  note  that  our  uni¬ 
form  grid  calculation;  arc  second-order  accurate.  Since  our  adaptive  calculations  arc  substantially 
more  accurate  for  a  given  computational  co.st,  the  i.ssuc  of  formal  truncation  error  is  moot. 
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7.  Data  Structures  and  Related  Issues 


From  the  outset  of  this  work  wc  adopted  a  retinement  strategy  that  would  allow  any  qua¬ 
drant  of  any  patch  to  be  refined  if  it  satisfied  some  criterion.  When  a  quadrant  is  refine,  each  mass 
control  volume  is  divided  into  four,  so  that  our  fundamental  data  objects  remain  identical,  even 
lliough  they  have  difi'erent  mesh  spacing  and  are  a.s.sociatcd  with  a  difTercnt  level  in  the  gnd 
hierarchy.  The  entire  patch  structure  that  we  generate  is  described  by  a  quad  tree. 


Multigrid  algorithms  for  nonlinear  PDFs  (full  approximation  storage  schemes)  calculate  a 
quantity  known  as  the  defect  as  an  integral  part  of  the  iteration  strategy.  One  definition  is 


xl’  =  //(' 


which  may  bo  thought  of  as  an  approximation  to  the  truncation  error  on  the  coarse  grid.  Other 
investigators  have  used  simpler  error  estimators  ba.sed  on  (.say)  velocity  gradients  or  curvature; 
however,  we  prefer  to  use  the  defect,  since  it  docs  not  bias  particular  terms  in  the  dillerenlial 
equation. 

Throughout  tliis  paper  we  have  refined  any  quadrant  tliat  satisfies 
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where  nVis  is  the  number  of  quadrants  that  arc  visible  (i.c.,  unrefined).  Normally,  we  take 
power  =  1.0.  However,  tiiis  criterion  is  not  crucial  to  our  algorithm,  and  other  criteria  can  be 
substituted  with  minimal  effort. 

In  addition  to  tfiis  criterion,  two  other  mcchani.sms  arc  u.scd  to  control  refinement.  To 
ensure  that  our  grid  ratio  docs  not  grow  unreasonably,  wc  force  a  refinement  if  die  ratio  is 
(locally)  greater  tlvan  two.  Moreover,  in  the  immediate  neighborhood  of  a  re-entrant  comer,  if  wc 
subdivide  one  quadrant  that  touehe.s  the  comer,  then  wc  refine  all  three  quadrants. 

A  number  of  data  structures  arc  required  to  support  the  various  numerical  processes  that 
have  to  be  performed  during  our  iterations.  The  proper  design  of  die.se  data  structures  enables 
reasonably  simple  implementation  of  each  of  the  numerical  operations,  even  though  the  arrange¬ 
ment  of  the  data  is  now  more  complex.  Typical  lists  of  patch  descriptors  include  the  following; 

•  The  NEXT  patch  at  level  N .  (Linking  the  patches  on  each  level  this  way  allows  all  opera¬ 
tions  to  be  done  by  li.st  traversal,  rather  than  by  a  recursive  tree  walk. ) 

•  The  XEIGUBCHiS  of  each  patch  (simplifying  data  exchange  during  relaxation  and  interpo¬ 
lation  operations' 


10 


•  The  identifiers  of  the  children  of  each  patch  (KIDS)  (necessary  for  prolongation). 

•  The  AUNTS  (neighbors  of  Ihe  parents)  of  each  patch  (required  for  restriction). 

Boundary  conditions  and  geometric  de.scripiors  arc  also  maintained  for  each  patch,  enabling 
us,  in  principle,  to  handle  curv'cd  boundaries  or  stretched  meshes  for  resolution  of  boundary 
layers. 

With  this  infonnation,  the  basic  numerical  operations  can  be  implemented  as  loops  going 
along  all  the  patches  in  a  particular  grid  level.  A  typical  code  fragment  could  be 

ID  =  LevPt(N’)  !  First  patch  on  Level  N 

DO  I  =  l.LevCnFN')  !  All  patches  on  Level  N 

NumjDp(Patch(lD].Patch( Kld(ID)))  !  Perform  numerical  process 

ID  =  N'ext(ID)  !  Next  patch  at  this  level 

FNDDO 

where  Num  (^p  is  one  of  snioothing,  prolongation,  restriction,  exchange  of  edge  data,  or  applica¬ 
tion  of  boundary  conditions.  It  is  clear  tiiat  the  normal  muliigrid  cycling  controls  determine  tfic 
invocation  of  these  code  fragments. 

The  natural  order  is  to  process  patches  in  the  order  of  their  creation  (which  determines  tlieir 
place  in  the  list).  All  processes  (except  smoothing)  arc  independent  of  the  ordering  of  the  points. 
Experiments  have  shown  that  the  smoothing  rate  of  our  relaxation  process  is  insensitive  to  the 
ordering. 

The  aspect  ratio  of  all  cells  remains  constant  throughout  the  refinement  process.  Of  course, 
grids  built  of  rectangles  can  easily  be  used  by  appropriate  definition  of  the  original  grid.  How¬ 
ever,  the  smoothing  rate  degrades  if  ihc  a.spccl  ratio  is  very  different  from  unity.  We  arc  consider¬ 
ing  using  line  relaxation  to  avoid  this  effect,  but  some  kind  of  patch  sorting  is  required  before  this 
strategy  can  be  implcmen'cd. 

S.  Numerical  Results 

We  w  ish  to  focus  attention  on  three  main  points.  First,  the  rate  of  convergence  and  accu¬ 
racy  that  we  achieve  on  our  automatically  adapted  grids  must  be  comparable  with  our  standard 
multigrid  calculations.  .Second,  our  local  truncation  error  estimate  xU  ,  which  we  u.sc  to  control 
our  local  rclinement  strategics,  must  give  us  reasonable  patterns.  Third,  the  accuracy  we  achieve 
mu,;t  be  comparable  with  that  obtained  by  other  benchmarks.  Additionally,  we  shall  discuss 
briefly  our  interface  treatment  and  dcmons  ralc  the  geometric  flexibility  of  the  algnrith.m.  How¬ 
ever.  we  have  not  attempted  to  develop  optimal  strategies  for  all  aspects  of  this  work,  but  rather  to 
validate  the  basic  components. 


8.1  The  Driven-Cavity  Problem 

The  first  example  we  consider  is  the  driven-cavity  problem.  This  widely  known  problem  is 
well  documented  in  the  literature.  We  have  chosen  this  problem  because  it  is  geometrically  sim¬ 
ple  and  the  solution  is  well  understood.  However,  the  presence  of  singularities  in  the  solution 
clearly  is  going  to  cause  any  adaptive  code  difficulties.  This  fact  should  be  borne  in  mind  when 
considering  these  results. 

The  problem  domain  is  a  unit  square,  and  the  flow  is  driven  by  specifying  u  =  1  and  v  =  0 
on  the  top  wall.  On  the  other  walls,  we  specify  u  =  v  =  0.  TTierc  arc  singularities  in  the  horizon¬ 
tal  velocity  in  the  top  comers  of  the  domain. 


Table  1  (a).  Effect  of  Different  Interface  Treatments  with  Re  =  1 , 
Fixed  Refinement  Patterns,  and  First-Order  Boundary  Treatment 
(timings  on  a  Sun/3  with  Weitek  floating-point  accelerator) 


Grid 

16 

16+5 
16+5+5 
16+5+5+5 
1 6+5 +5 +5 +5 
512x512 


Table  1  (b).  Effect  of  Different  Interface  Treatments  with  Re  =  400, 
Fixed  Refinement  Patterns,  and  First-Order  Boundary  Treatment 


Linear  Interface 

Quadratic  Interface 

Grid 

Tk. 

Error 

Th. 

Error 

# 

64  X  64 

2.549(-2) 

5.2(-4) 

2.549(-2) 

5.2(-4) 

■a 

64  X  64+5 

2.503(-2) 

0.6(-4) 

12 

2.493(-2) 

■9 

64  X  64+10 

2.503(-2) 

0.6(-4) 

11 

2.500(-2) 

0.3(-4) 

12 

128 X  128 

2.479(-2) 

1.8(-4) 

11 

2.479(-2) 

1.8(-4) 

11 

128  X  128+5 

2.491  (-2) 

0.6(-4) 

10 

2.479(-2) 

1.8(-4) 

11 

128  X  128+10 

2.491  (-2) 

0.6(-4) 

2.487(-2) 

10 

64  X  64+5+5 

2.5l5(-2) 

1.8(-4) 

2.493(-2) 

0.4(-4) 

11 

Notes;  64  x  64  and  128  x  128  denote  uniform-grid  calculations. 

64  X  64+5  has  the  five  rows  of  mass  control  volumes  at  the  top  of  the  domain  subdivided 
(each  cell  becoming  four  finer  cells). 

64  X  64+5+5  has  the  top  five  rows  of  cells  of  the  64  x  64+5  grid  subdivided. 

Other  definitions  are  obvious  generalizations  of  this. 


First,  wc  consider  ihc  effect  of  the  interface  treatment  on  the  accuracy  of  the  solution. 
Tables  1(a)  and  1(b)  show  the  wall  shear  stress  and  the  associated  errors  at  the  midpoint  of  the 
wall.  The  wall  shear  value  shown  in  the  tables  was  calculated  by  using  third-order,  one-sided 
interpolation.  The  error  was  measured  by  taking  the  value  calculated  on  a  512  x  512  grid  as 
exact.  (The  discretized  equations  used  a  one-sided,  first-order  approximation  for  the  wall  shear, 
for  simplicity.)  For  the  purposes  of  this  comparison  only,  we  have  forced  the  refinement  patterns. 
The  grid  refinements,  explained  in  the  table  notes,  were  chosen  to  allow  reasonable  representa¬ 
tions  of  (he  boundary  layer. 

From  these  tests  (and  others  at  even  higher  Reynolds  numbers),  we  deduced  that  second- 
order  interpolation  was  adequate  for  our  fully  adaptive  strategy.  (As  we  expected,  the  quadratic 
formulae  arc  somewhat  more  accurate,  but  the  gains  do  not  seem  to  justify  the  additional  com¬ 
plexity.)  Clearly,  this  approach  differs  from  that  described  in,  for  example,  reference  [6].  Our 
goal  is  to  introduce  refinements  that  efficiently  reduce  the  errors  associated  with  the  existing  grid 
calculations,  rather  than  to  maintain  a  particular  order  of  truncation  error.  Formal  truncation  error 
analysis  and  its  relationship  to  actual  errors  is  a  complex  issue  in  this  context.  We  are  not  claim¬ 
ing  that  our  adaptive  di.scrctization  is  second  order.  However,  our  uniform  grid  calculations  are, 
in  general,  second  order.  Throughout  this  section  wc  compare  the  adaptive  and  the  uniform  (i.c., 
second-order)  .strategies  on  the  basis  of  speed,  accuracy,  and  .storage. 


Table  1(c)  shows  the  performance  of  the  fully  automatic  adaptive  code,  relative  to  uniform 
grid  calculations.  The  Reynolds  number  is  400.  The  errors  have  be  m  estimated  by  considering 
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the  wall  shear  stresses  at  the  midpoint  of  the  moving  lid.  We  have,  of  course,  deliberately  chosen 
a  measure  that  is  not  in  a  neighborhood  of  the  comer  singularities.  Winters  and  Cliffe  [5]  have 
used  two  independent  formulations  for  their  computations  (both  primitive  variable  and  stream- 
function/vorticity).  The  values  they  give,  Ix=.5  =  2.661  (-2)  or  2.481  (-2),  can  be  used  as  an  esti¬ 
mate  of  the  accuracy  of  their  figures.  It  is  clear  that  both  our  uniform  and  adapted  grids  converge 
to  similar  answers.  However,  our  present  refinement  strategy  concentrates  on  the  singularities  in 
the  comers.  The  wall  shear  in  the  middle  of  the  lid  does  not  app)ear  to  be  especially  sensitive  to 
the  comer  treatment.  (This  seems  reasonable:  it  is  in  the  middle  of  a  boundary  layer  driven  by  the 
moving  wall.)  Therefore,  the  adaptive  calculations  will  be  somewhat  less  effective  than  uniform 
ones  for  this  particular  measure  of  accuracy  (see  Table  1(c)). 


Table  1(c).  Driven-Cavity  Problem:  Wall  Shear  Stress  at  Midpoint  of  Moving  Lid  with  Re  =  400 


Minimum 

Grid  Size 

Level 

Uniform  Calc.  (2nd  order) 
tw  li=.5  Time  (sec)  #  Patch 

Adaptive  Calc.  (2nd  order) 

Th- ljt=.5  Time  (sec)  #  Patch 

1/16 

3 

0.03601 

6 

21 

0.03601 

6 

21 

1/32 

4 

0.02825 

17 

85 

0.02937 

17 

43 

1/64 

5 

0.02517 

46 

341 

0.02793 

34 

101 

1/128 

6 

0.02505 

109 

1365 

0.02633 

66 

220 

1/256 

7 

0.02508 

828 

5461 

0.02554 

126 

462 

1/512 

8 

- 

- 

- 

0.02539 

277 

954 

1/256^ 

7 

- 

- 

- 

0.02554 

85 

462 

1/512<> 

8 

- 

- 

- 

0.02539 

213 

954 

1/256* 

7 

- 

- 

- 

0.02541 

99 

540 

1/512* 

8 

- 

- 

- 

0.02522 

256 

1229 

l/256<^ 

7 

- 

- 

- 

0.02516 

293 

1637 

Notes: 

F(2,2)  cycles,  tol=1.0(-7),rcf-power=0.35. 

Timings  on  an  IBM  Rise  Systcm/6000®  (model  530). 

"  results  obtained  with  the  ad  hoc  cycling  strategy  described  at  the  end  of  Section  8.2. 

results  obtained  with  ref-power  =  0.2  and  the  ad  hoc  cycling  strategy  (ref-power  =  0.35 
otherwise). 

‘^results  obtained  with  the  ad  hoc  cycling  strategy  described  below  on  an  Adapt(5,7)  grid. 
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The  vertical  velocity  extrema  at  the  horizontal  mid-plane  are  an  alternative  error  measure. 
In  Table  1(d)  we  compare  adaptive  and  uniform  calculations  on  the  basis  of  these  quantities. 
Again,  we  see  the  importance  of  the  refinement  strategy.  However,  the  Adapt(5,7)  calculation  has 
accuracy  comparable  with  the  uniform  seven-level  calculation  (at  least  for  the  minimum  vertical 
velocity)  with  considerable  savings  in  epu  time  and  computer  storage.  Physically,  one  could 
argue  that  the  downturn  in  the  horizontal  boundary  layer  is  caused  by  the  pressure  maximum  at 
the  top  right-hand  comer  of  the  domain.  Our  adaptive  strategy  has  captured  this  downturn  and 
has  also  refined  the  mesh  downwards,  following  this  “jet.”  In  Figure  5,  we  show  a  seven-level 
adaptive  grid  and  the  resulting  flow  visualization  consisting  of  the  level  curves  for  the  stream 
function. 


Table  1(d).  Driven-Cavity  Problem;  Velocity  Extrema  at  the  Horizontal  Midplane 


Minimum 
Grid  Size 

Level 

Comment 

Midplane  Y=0.5 

^  min  T  max 

Time  (sec) 

#  Patch 

3 

Uniform 

-0.2931 

6 

21 

4 

Uniform 

-0.3923 

KB 

17 

85 

1/64 

5 

Uniform 

-0.4396 

0.2946 

46 

341 

1/128 

6 

Uniform 

-0.4518 

0.3021 

109 

1365 

1/256 

Uniform 

-0.4534 

0.3033 

828 

5461 

1/25  6<> 

BH 

Adapt(3,7) 

-0.4388 

0.2852 

85 

462 

1/2560 

Bi 

Adapt(5,7) 

-0.4523 

0.3020 

293 

1637 

Notes; 

Validation  and  serial  performance  for  the  adaptive  code;  vertical  velocity  extrema  across  the  hor¬ 
izontal  midplane. 

Adapt(n,m)  defines  a  calculation  with  n  levels  of  uniform  refinement  and  then  adaptive 
refinement  up  to  level  m .  (In  this  table  we  had  ref-power  =  0.35.) 

Timings  on  an  IBM  Rise  Systcm/6000®  (model  530). 

results  obtained  with  the  ad  hoc  cycling  strategy  described  at  the  end  of  Section  8.2. 


Table  1(e)  shows  the  rates  of  convergence  for  the  driven-cavity  problem  on  uniform  and 
adapted  grids.  As  is  typical  of  multigrid  algorithms  on  nonlinear  problems,  the  convergence  rate 
varies  somewhat  as  the  number  of  levels  increases  but  remains  bounded  weU  away  from  one  uni¬ 
formly  in  the  number  of  grid  levels. 

Of  greater  significance,  we  note  that  the  convergence  rates  achieved  on  the  adapted  meshes 
are  comparable  to  the  uniform  grid  rates.  One  would  not  necessarily  expect  this  result  on  theoret¬ 
ical  grounds;  it  amounts  to  a  pleasant  surprise.  The  calculations  here  are  terminated  on  each  grid 
when  the  norm  of  the  residual  is  less  than  10.0(-07).  The  reduction  factors  given  in  the  table  arc 
the  final  reduction  factors',  in  some  cases  they  have  not  reached  their  asymptotic  values. 
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Each  of  our  test  cases  contains  some  form  of  singularity  (or  singularities).  We  are  not  treat¬ 
ing  these  in  any  special  fashion;  they  are  automatically  handled  by  our  refinement  strategy.  In 
general,  for  complex  flows,  it  is  impossible  to  predict  all  possible  singularities;  therefore, 
methods  that  rely  on  accurate  and  special  treatment  of  these  regions  are  not  sufficiently  flexible. 
Our  solutions  are  obviously  going  to  be  poor  in  small  neighborhoods  of  the  singularities,  but  this 
effect  is  locai,  moreover,  careful  mesh  grading  allows  a  solution  to  be  obtained  without  degrada¬ 
tion  of  the  order  of  accuracy  (see,  for  example,  (3)  and  [8]). 


Table  1(e).  Driven-Cavity  Problem: 

Rate  of  Convergence  on  Adaptive  and  Nonadaptive  Grids,  with  Re  =  400 


Grid 

Level 

Uniform  Calculation  (2nd) 

Total  Fine  No. 

Pats.  Pats.  Cycles  pf 

Adapted  Calculation  (2nd) 

Total  Fine  No. 

Pats.  Pats.  Cycles  pf 

5 

4 

23 

.51 

5 

4 

23 

.51 

21 

16 

20 

.47 

21 

16 

20 

.47 

85 

64 

24 

.52 

43 

22 

19 

.42 

341 

256 

22 

.55 

101 

52 

17 

.23 

1365 

1024 

15 

.31 

220 

102 

16 

.47 

5461 

4096 

20 

.62 

462 

169 

15 

.42 

8 

- 

- 

" 

954 

321 

17 

.52 

Notes: 

F(2,2)  cycles,  tol.  =  1.0(-7),  ref-power  =  0.35. 
Second-order  boundary  conditions. 


8.2  The  Flow  over  a  Backward-Facing  Step  Problem 

The  flow  over  a  backward-facing  step  is  another  geometrically  simple  and  well-documented 
problem.  As  is  usual,  we  have  specified  a  parabolic  inlet  profile  (fully  developed  flow),  and  we 
have  fully  developed  outlet  conditions.  No  slip  conditions  are  imposed  at  all  walls.  In  this 
configuration,  we  calculate  the  flow  upstream  of  the  step,  thereby  solving  for  the  flow  around  the 
step  and  not  ignoring  the  singularity  at  the  comer.  In  contrast  to  [15],  important  features  of  this 
flow  are  the  recirculation  zone  behind  the  step  and  the  velocity  profiles  downstream  of  the  step. 

In  Table  2(a)  we  give  the  convergence  rates  for  the  uniform-grid  and  adapted-grid  calcula¬ 
tions.  Little  degradation  is  observed  in  the  adapted-grid  computations,  despite  the  complexity  of 
the  problem  domain.  However,  the  uniform-grid  calculations  at  level  3  converge  more  quickly 
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because  the  diffusion  terms  become  more  important  than  the  convective  terms  as  the  mesh-size 
decreases.  This  effect  is  not  seen  on  the  adaptive  calculations  because  some  parts  of  the  coarser 
grid  remain  visible. 


Table  2(a).  Rate  of  Convergence  for  the  Backward-Facing  Step  Problem: 
Uniform  and  Adaptive  Calculations 


Total 

Pats. 

Uniform  Calculation 

Fine  No. 

Pats.  Cycles 

Total 

Pats. 

Adapted  Calculation 
Fine  No. 

Pats.  Cycles 

MM 

144 

16 

.52 

180 

144 

16 

.52 

576 

9 

.24 

240 

60 

11 

.41 

H 

n 

- 

- 

- 

391 

114 

10 

.40 

hh 

■ 

— 

— 

— 

525 

19 

14 

.42 

■i 

“ 

- 

- 

544 

19 

15 

.42 

Notes:  F(2,2)  cycles  with  Re  =  150,  power  =  2.0,  and  tol.  =  10.0(-6). 


Table  2(b).  Results  for  the  Backward-Facing  Step: 
Uniform  and  Adaptive  Calculations  and  Benchmark 


Velocity 

- 

Cliffe 

(2,2) 

(3.3) 

(2.4) 

(2,10) 

U  max 

1.6  step-ht 

0.904 

0.906 

0.905 

4.0  step-ht 

0.725 

0.722 

0.726 

0.725 

f/ min 

1 .6  Step-ht 

-0.102 

-0.103 

-0.103 

-0.102 

-0.101 

4.0  step-ht 

-0.049 

-0.036 

-0.047 

-0.039 

-0.038 

Umax 

1 .6  step-ht 

0.010 

0.010 

0.010 

0.010 

0.010 

4.0  step-ht 

0.0 

0.0 

0.0 

0.0 

0.0 

U  min 

1.6  step-ht 

-0.070 

-0.068 

-0.070 

-0.071 

-0.071 

4.0  step-ht 

-0.091 

-0.077 

-0.087 

-0.089 

-0.089 

No.  patches 

— 

180 

756 

391 

544 

Time  (sec) 

- 

73 

77 

113 

- 

Total  Time  (sec) 

73 

77 

278 

Notes:  (2,2):  uniform  2-level  calculation  with  mesh  size  =  step-height/8. 

(3,3):  uniform  3-level  calculation  with  mesh  size  =  step-height/16. 

(2,N):  adapted  calculation  with  maximum  mesh  size  =  step-height/8  and  minimum  mesh  size 
=  step-height/2(^+*\ 


To  assess  the  accuracy  of  our  approach,  we  have  included  uniform,  adapted,  and  benchmark 
solutions;  see  Table  2(b).  The  benchmark  solutions  were  taken  from  reference  [14].  Clearly,  all 
the  sets  of  calculations  are  in  agreement.  Moreover,  the  adaptive  calculations  with  two  levels  of 
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automatic  refinement  (shown  as  (2,4)  in  the  table)  use  only  half  the  storage  of  the  uniform  3-level 
calculations  for  comparable  accuracy.  A  weak  singularity  in  the  neighborhood  of  the  comer  inhi¬ 
bits  further  improvement  of  the  solution  with  our  current  refinement  strategy.  However,  there  is 
no  degradation  in  the  solution  if  we  allow  more  refinement  (we  allowed  a  further  six  levels). 
Therefore,  our  interface  treatment  is  reasonable. 

We  show  two  sets  of  flow  visualizations  and  refinement  patterns  for  this  configuration.  The 
first  set  (Figure  6(a))  was  obtained  by  using  a  naive  treatment  of  the  inlet  and  first-order  treat¬ 
ments  of  the  walls.  The  code  has  automatically  detected  this  and  refined  the  pattern  in  the  areas 
where  the  approximation  is  poor.  When  we  improved  these  two  aspects  so  that  we  had  a  con¬ 
sistently  second-order  treatment  throughout  the  domain,  we  obtained  the  refinement  patterns 
shown  in  Figure  6(b)  (the  values  used  in  Tabic  2(b)  used  the  improved  treatments). 

Our  coarsest-grid  treatment  simply  consists  of  many  sweeps.  In  many  circumstances  this  is 
adequate.  However,  in  cases  where  there  are  many  patches  on  the  first  grid,  it  is  obvious  that 
various  forms  of  acceleration  could  be  used.  At  present,  these  have  not  been  implemented  since 
the  benefits  are  probably  marginal.  Moreover,  we  can  decrease  this  effect  by  making  minor 
modifications  to  the  multigrid  cycling  strategy. 

Another  factor  that  affects  our  adaptive  timings  is  that,  at  present,  we  solve  a  sequence  of 
flow  problems  so  that  we  have  reliable  error  estimates  that  are  “uncontaminated”  by  insufficient 
numerical  convergence;  these  are  the  basis  of  our  refinement  strategy.  It  may  be  possible  to  avoid 
this  .^cqucntial  procedure  by  doing  the  refinement  and  solution  concurrently  (in  some  sense). 
However,  we  have  avoided  this  strategy  in  the  present  study  since  the  convergence  errors  confuse 
the  refinement  patterns.  The  only  attempt  we  have  made  to  address  this  issue  is  the  ad  hoc 
strategy  referred  to  in  the  tables.  Here  we  do  the  following:  when  performing  an  Adapt {n  jn)  cal¬ 
culation  (uniform  refinement  on  n  grids,  then  adaptive  refinement  up  to  level  m),  we  reset  the 
cycling  tolerance  to  its  square  root  until  we  have  completed  the  preliminary  problems  and  we 
have  an  adapted  m -level  grid.  We  then  solve  the  m -level  problem  to  the  full  tolerance.  This 
strategy  decreases  the  set-up  time,  without  significantly  changing  the  final  adapted  grid  or  the 
computed  solution.  For  two  of  the  geometries  in  this  paper,  this  strategy  is  quite  effective.  How¬ 
ever,  for  the  backward-facing  step  geometry  we  were  able  to  reduce  the  time  for  the  Adapt (2A) 
calculation  only  to  180  seconds.  The  coarsest-grid  effects  still  dominate  when  we  have  such  a 
small  number  of  levels. 


8.3  The  Sudden  ExpansionI Contraction  Geometry 

The  third  example  is  new.  It  consists  of  a  sudden  expansion/contractiuii  ouu  was  cnostn  to 
show  the  flexibility  of  our  approach.  (The  choice  was  also  tempered  by  the  desire  to  have  an 
intuitive  feeling  for  the  correctness  of  our  results.)  The  geometry  is  shown  in  Figure  7.  At  the 
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extreme  left,  we  specify  a  parabolic  inlet  profile  for  the  horizontal  velocity  (and  v  =  01.  On  the 
right  we  specify  fully  developed  flow  conditions:  v  =  0  and  =  0.  Everywhere  else  we  have 

a  =  V  =  0.  This  flow  has  four  re-entrant  comers.  An  interesting  quantity  is  the  flux  across  the 
“entrance”  to  the  expansions  in  the  domain. 

The  rates  of  convergence  (Table  3)  agree  with  the  results  of  the  preceding  sections;  adap¬ 
tivity  does  not  slow  down  convergence.  Again  we  see  the  uniform-grid  calculations  improving  as 
the  mesh  becomes  finer,  because  the  elliptic  terms  become  dominant;  but  this  effect  is  not  seen  in 
the  adaptive  calculations,  because  sections  of  the  coarse  grids  remain  “visible”  and  limit  the  rate 
of  convergence.  The  flow  visualization  (Figure  7)  shows  that  we  are  obtaining  qualitatively  rea¬ 
sonable  results.  There  are  two  shear-driven  recirculation  cells  formed  behind  the  steps.  The 
streamlines  seem  to  be  quite  reasonable. 

The  streamlines  shown  in  the  figures  were  calculated  by  integrating  over  the  adapted  grid. 
The  use  of  this  adapted  mesh  was  necessary  because  the  conservation  equation  is  satisfied  here. 
',Ve  used  Cuussian  quadrature,  modified  to  capture  the  compound  grid  automatically.  Techniques 
that  ignored  this  fact  gave  disappointing  results. 

One  quantity  of  interest  is  the  transport  from  the  main  channel  into  the  two  expansions.  We 
have  measured  this  by  looking  at  the  vertical  velocity  extrema  along  the  “continuation”  line 
across  the  top  of  the  channel.  Both  the  uniform  and  adaptive  calculations  converge  to  the  same 
values;  however,  when  the  ad  hoc  cycling  strategy  is  used,  the  adaptive  calculations  are  faster  and 
use  less  memory. 
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Table  3.  Sudden  Expansion/Contraction  Geometry; 
Rate  of  Convergence  and  Vertical  Velocity  Extrema 


Unifonn  Calculation 

Adaptive  Calculation 

Minimum 

Grid  Size 

Vert.  Velocity- 

Max  Min 

Time 

(sec) 

No. 

Patches 

No. 

Cycles 

Vert.  Velocity 

Max  Min 

Time  (sec) 

Item  Cumul. 

No. 

Patches 

No. 

Cycles 

1/8 

0.0212 

-0.073 

6 

25 

11 

0.0212 

-0.073 

6 

6 

25 

11 

1/16 

0.0224 

-.088 

9 

105 

9 

0.0230 

-0.090 

11 

16 

55 

14 

1/32 

0.0233 

-.092 

22 

425 

8 

0.0243 

-0.094 

18 

34 

Ill 

12 

1/64 

0.0254 

-.099 

72 

1705 

8 

0.0254 

-0.100 

27 

61 

185 

14 

1/128 

... 

... 

- 

.. 

- 

0.0262 

-0.100 

44 

106 

283 

14 

l/64‘> 

... 

... 

- 

- 

- 

0.0257 

-0.100 

31 

46 

186 

14 

1/1280 

--- 

— 

0.0262 

-0.100 

44 

69 

284 

14 

Notes: 

Re  =  50,  F(2,2)  cycles,  tol=1.0(-7),  ref-power^  1.0. 
Timings  on  an  IBM  Rise  System/6000®  (model  530). 

Results  obtained  with  the  ad  hoc  cycling  .strategy. 


9.  Discussion 

There  are  three  reasons  why  one  is  interested  in  automatic  adaptive  gridding  strategies: 
speed,  storage,  and  reliability  of  the  solution.  The  first  item  has  been  dealt  with  earlier. 

The  issue  of  storage  is  fairly  simple,  and  we  have  been  able  to  achieve  significant  savings  in 
all  cases  of  interest.  Here  we  simply  add  the  comment  that,  in  three-dimensional  problems,  the 
question  of  efficient  storage  utilization  is  significantly  more  important  and  the  benefits  from  adap¬ 
tive  gridding  correspondingly  greater. 

The  issue  of  reliability  of  the  computed  solution  is  complex.  Our  present  refinement 
strategy  is  somewhat  ad  hoc,  and  a  detailed  study  of  refinement  strategies  is  beyond  the  scope  of 
this  paper.  However,  our  examples  illustrate  three  types  of  error:  (1)  those  induced  by  some  type 
of  singularity,  (2)  those  caused  by  comparatively  poor  numerical  treatments  (e.g.,  low-order 
boundary  formulae),  and  (3)  those  caused  by  insufficient  grid  resolution.  The  best  that  one  can 
hope  to  do  in  the  presence  of  singularities  is  to  treat  them  sufficiently  well  that  their  effect  is  local 
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(unless  one  uses  some  semianalylic  approach).  Our  code  does  this.  Moreover,  the  other  two 
effects  are  detected  and  automatically  reduced  by  our  refinement  strategy  until  the  errors  are  dom¬ 
inated  by  singularity  effccLs. 


The  question  of  how  a  numerical  solution  varies  with  respect  to  “small”  changes  in  the  grid 
is  not  well  understood,  nor  is  the  global  effect  of  singularities  in  the  domain.  Attention  must  be 
given  to  these  points  before  we  can  achieve  “full  reliability.”  However,  it  is  interesting  to  note 
that  the  problems  with  singularities  arc  not  restricted  to  error  estimators  based  on  the  defect  (or 
the  local  truncation  error).  Richardson  extrapolation  techniques  and  simpler  gradient  estimators 
can  have  the  same  problems. 


Throughout  this  work  we  have  measured  errors  by  estimating  a  continuous  solution  and 
defining  errors  with  respect  to  this.  We  have  the  following  equation: 


L* 


where  is  the  discrete  error.  Since  the  equation  is  linear,  it  can  be  solved  with  comparatively 
little  extra  work,  given  that  the  operator  and  the  gridding  are  already  defined.  This  approach 
offers  possibilities  of  improving  our  refinement  strategies.  However,  may  contain  nonlocal 
errors,  thereby  complicating  its  use  in  a  refinement  strategy. 


10.  Conclusions 

As  is  customary  with  this  type  of  work,  we  have  demonstrated  the  effectiveness  and  accu¬ 
racy  of  our  techniques  by  considering  a  few,  well-known  model  problems. 

We  have  shown  that  the  rate  of  convergence  we  ob.serve  for  the  adaptive  calculations  does 
not  degrade  significantly  compared  to  the  uniform  ones.  Our  estimate  of  the  local  truncation  error 
(x/!0  gives  us  qualitatively  correct  refinement  patterns,  which  are  capable  of  detecting  deficiencies 
in  the  numerical  treatment  as  well  as  in  the  grid.  The  accuracy  of  our  results  is  comparable  with 
that  of  other  benchmarks  for  both  the  uniform  and  the  adaptive  calculations.  Our  current  adaptive 
strategy  will  deal  with  all  three  clas.ses  of  error.  However,  it  continues  to  “home  in”  on  singulari¬ 
ties,  even  though  this  is  not  always  an  optimal  strategy. 

Additionally,  we  have  shown  that  the  interpolation  that  we  use  at  interfaces  is  good  and  that 
we  are  easily  able  to  handle  a  range  of  geometries.  For  some  of  the  cases  we  have  considered,  we 
have  been  able  to  achieve  significant  gains  in  speed  and  storage  compared  to  uniform  grid  calcu¬ 
lations.  However,  we  have  shown  that  this  is  not  always  the  case. 


11.  Future  Work 


In  liic  current  a  uk  wc  showed  that  we  can  separate  the  three  components-the  multigrid 
algoritiim,  the  adaptive  strategy,  and  the  numencal  approximations-as  was  our  design  goal.  Now 
tliat  we  have  done  Uiis,  we  can  modify  and  improve  each  of  Uicse  aspects  separately.  In  particu¬ 
lar,  wc  are  looking  at  eflicient  ways  to  improve  the  generality,  .scope,  and  accuracy  of  Ihe  algo- 
nlli.m.  We  are  afso  focu-ing  our  attention  on  parallcli/aiion  i.ssucs.  A  shared -memory'  version  of 
this  code  is  ainmug,  ai'd  we  are  atout  to  look  at  implementation  on  scalable  distributed-memory' 
architecture^-,  such  as  Uie  Intcl0  Touchstone  DELTA  system. 
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Figure  1(a).  Driven  Cavity: 
Solution  Domain  and  Boundary 
Conditions 


Figure  1(b).  Staggered  Grid:  Location 
of  Variables.  Bracketed  variables  are 
associated  with  adjacent  cells. 
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Figure  3(a).  The  Canonic  Problem:  One  Coarse-Grid  Patch  Interfacing 
with  Four  Fine-Grid  Patches.  This  tiling  is  sufficient  to  generate 
all  configurations  of  interest.  We  show  the  location  of  "genuine" 
unknowns  on  the  compound  grid.  At  these  locations  we  solve  approxi¬ 
mations  of  the  differential  equations.  (Coarse-grid  variables  are  in 
upper  case,  fine-grid  ones  in  lower  case.) 
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Figure  3(b).  The  Canonic  Problem  Modified  to  Show  the  Dummy  Unknowns 
Introduced  for  Computational  Convenience.  The  variables  shown  here 
are  not  unknowns  in  our  system  of  equations  but  dummy  variables  that 
are  intioduced  for  convenience.  Special,  additional  equations  are 
introduced  for  them.  (Coarse-grid  variables  are  in  upper  case,  fine- 
grid  ones  in  lower  case.) 
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Figure  4.  Prolongation  at  Interfaces.  (Coarse-grid 
unknowns  are  shown  in  upper  case,  fine-grid  ones  in 
lower  case.) 


Figure  5.  The  Driven  Cavity  Problem;  Flow  visualization,  an  adapted  grid  and  problem 
definition.  On  the  top  wall  we  drive  the  fiow  by  specifying  u  =  1  (and  v  =  0).  On  the  other 
walls  we  have  m  =  v  =  0.  There  are  singularities  in  the  horizontal  velocity  in  the  two  top 
comers  of  the  domain. 


(h) 


Figure  6.  The  Backward-Facing  Step  Problem:  Flow  visualization,  an  adapted  grid  and 
problem  definition. 

(a)  First-order  boundary  conditions  and  standard  inlet  treatment, 

(b)  Second-order  boundary  conditions  and  modified  inlet  treatment. 

At  the  extreme  left  we  specify  a  parabolic  inlet  profile  for  the  horizontal  velocity  (and 
V  =  0).  On  the  right  we  specify  fully  developed  flow  conditions:  v  =  0  and  =  0.  On  all 

other  boundaries  we  have  m  =  v  =  0.  Notice  that  we  are  solving  the  flow  around  the  step 
and  not  ignoring  the  singularity  at  the  comer.  This  is  a  significantly  more  difficult  problem 
than  considering  a  rectangular  domain  only. 
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Figure  7.  The  Cross  Geometry;  Flow  visualization,  an  adapted  grid,  and  problem  definition 
At  the  extreme  left  we  specify  a  parabolic  inlet  profile  for  the  horizontal  velocity  (and  v  =  0).  On 
the  right  we  specify  fully  developed  flow  conditions;  v  =  0  and  ||-  =  0.  On  all  oih.  r  boundaries 

we  have  u  =  V  =  0.  This  flow  has  four  re-entrant  comers.  An  interesting  quantity  is  the  flux 
across  the  entrance”  to  the  expansions  in  the  domain.  Peak  vclociUcs  arc  shown  in  Table  3. 


IVUNSA 


Report  Documentation  Page 


1  Report  No  2.  Government  Accession  No 

NASA  CR- 187595 
ICASE  Report  No.  91-36 
4  Title  and  Subtitle 

A  DYNAMICALLY  ADAPTIVE  MULTIGRID  ALGORITHM  FOR  THE 
INCOMPRESSIBLE  NAVIER-STOKES  EQUATIONS  —  VALIDATION 
AND  MODEL  PROBLEMS 


3  Recipient's  Catalog  No 

5  Report  Date 

June  1991 

6  Performing  Organization  Code 


7  Authorls) 

C.  P.  Thompson 
G.  K.  Leaf 
J.  Van  Rosendale 

9.  Performing  Organization  Name  and  Address 

Institute  for  Computer  Applications  in  Science 
and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 
Hampton,  VA  23665-5225 


12  Sponsoring  Agency  Name  and  Address 
National  Aeronautics  and  Space  Administration 
Langley  Research  Center 
Hampton,  VA  23665-5225 


IS.  Supplementary  Notes 
Langley  Technical  Monitor; 
Michael  F.  Card 


8.  Performing  Organization  Repon  No 

I  91-36 

10  Work  Unit  No. 

I  505-90-52-01 

[  11  Contract  or  Grant  No 

NASl-18605 

13.  Type  of  Report  and  Period  Covered 
Contractor  Report 


14  Sponsoring  AQency  Code 


Submitted  to  J,  Applied  Numerical 
Mathematics 


Final  Report 


16  Abstract 

We  describe  an  algorithm  for  the  solution  of  the  laminar,  incompressible  Navier-Stokes  equations.  The  basic 
algorithm  is  a  multigrid  method  based  on  a  robust,  box-based  smoothing  step.  Its  most  important  feature  is  the 
incorporation  of  automatic,  dynamic  mesh  refinement.  Using  an  approximation  to  the  local  truncation  error  to 
control  the  refinement,  we  use  a  form  of  domain  decomposition  to  introduce  patches  of  finer  grid  wherever  they 
[  are  needed  to  ensure  an  accurate  solution.  This  refinement  strategy  is  completely  local:  regions  that  satisfy  our 
tolerance  are  unmodified,  except  when  they  must  be  refined  to  maintain  reasonable  mesh  ratios  This  locality  has 
the  important  consequence  that  boundary  layers  and  other  regions  of  sharp  transition  do  not  “steal”  mesh  points 
from  surrounding  regions  of  smooth  flow,  in  contrail  to  moving  mesh  strategies  where  such  “stealing"  is  inevitable. 

Our  algorithm  supports  generalized  simple  domains,  that  is,  any  domain  defined  by  horizontal  and  vertical  lines. 
This  generality  is  a  natural  consequence  of  our  domain  decomposition  approach.  We  base  our  program  on  a  standard 
staggered-grid  formulation  of  the  Navier-Stokes  equations  for  robustness  and  efficiency.  To  ensure  discrete  mass 
conservation,  we  have  introduced  special  grid  transfer  operators  at  grid  interfaces  in  the  multigrid  algorithm.  While 
these  operators  complicate  the  algorithm  somewhat,  our  approach  results  in  exact  mass  conservation  and  rapid 
convergence. 

In  this  paper,  we  present  results  for  three  model  problems:  the  driven-cavity,  a  backward-facing  step,  and 
a  sudden  expansion/contraction.  Our  approach  obtains  convergence  rates  that  are  independent  of  grid  size  and 
that  compare  favorably  with  other  multigrid  algorithms.  We  also  compare  the  accuracy  of  our  results  with  other 
benchmark  results. 
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